AFRL-RV-PS- 

TR-2015-0182 


AFRL-RV-PS- 

TR-2015-0182 


CHARACTERIZATION  OF  N ON -LINEARIZED 
SPACECRAFT  RELATIVE  MOTION  USING 
NONLINEAR  NORMAL  MODES 


Eric  A.  Butcher 


University  of  Arizona 
888  N  Euclid  Ave 
Tucson,  AZ  85719-4824 


20  Apr  2016 


Final  Report 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  IS  UNLIMITED. 


AIR  FORCE  RESEARCH  LABORATORY 

Space  Vehicles  Directorate 

3550  Aberdeen  Ave  SE 

AIR  FORCE  MATERIEL  COMMAND 

KIRTLAND  AIR  FORCE  BASE,  NM  87117-5776 


DTIC  COPY 

NOTICE  AND  SIGNATURE  PAGE 


Using  Government  drawings,  specifications,  or  other  data  included  in  this  document  for 
any  purpose  other  than  Government  procurement  does  not  in  any  way  obligate  the  U.S. 
Government.  The  fact  that  the  Government  formulated  or  supplied  the  drawings, 
specifications,  or  other  data  does  not  license  the  holder  or  any  other  person  or  corporation; 
or  convey  any  rights  or  permission  to  manufacture,  use,  or  sell  any  patented  invention  that 
may  relate  to  them. 


This  report  is  the  result  of  contracted  fundamental  research  deemed  exempt  from  public  affairs 
security  and  policy  review  in  accordance  with  SAF/AQR  memorandum  dated  10  Dec  08  and 
AFRL/CA  policy  clarification  memorandum  dated  16  Jan  09.  This  report  is  available  to  the  general 
public,  including  foreign  nationals.  Copies  may  be  obtained  from  the  Defense  Technical 
Information  Center  (DTIC)  (http://www.dtic.mil). 


AFRL-RV-PS-TR-20 15-0 182  HAS  BEEN  REVIEWED  AND  IS  APPROVED  FOR 
PUBLICATION  IN  ACCORDANCE  WITH  ASSIGNED  DISTRIBUTION  STATEMENT. 


//SIGNED// 

THOMAS  LOVELL 
Program  Manager 


//SIGNED// 

PAUL  HAUSGEN,  Ph.D. 

Technical  Advisor,  Spacecraft  Component  Technology 


//SIGNED// 

JOHN  BEAUCHEMIN 

Chief  Engineer,  Spacecraft  Technology  Division 
Space  Vehicles  Directorate 


This  report  is  published  in  the  interest  of  scientific  and  technical  information  exchange,  and  its 
publication  does  not  constitute  the  Government’s  approval  or  disapproval  of  its  ideas  or  findings. 


Approved  for  public  release;  distribution  is  unlimited. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and  maintaining  the 
data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing 
this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202- 
4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently 
valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  (DD-MM-YY)  2.  REPORT  TYPE  3.  DATES  COVERED  (From  -  To) 

20-04-20 1 6  Final  Report  3  Sep  2014-03  Mar  2016 


4.  TITLE  AND  SUBTITLE  5a.  CONTRACT  NUMBER 

Characterization  of  Non-Linearized  Spacecraft  Relative  Motion  using 

Nonlinear  Normal  Modes  FA9453-14-1-0350 

5b.  GRANT  NUMBER 


6.  AUTHOR(S) 

Eric  A.  Butcher 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Arizona 
888  N  Euclid  Ave 
Tucson,  AZ  85719-4824 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory 
Space  Vehicles  Directorate 
3550  Aberdeen  Ave.  SE 
Kirtland  AFB,  NM  87117-5776 


12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 


5c.  PROGRAM  ELEMENT  NUMBER 

62601F 


5d.  PROJECT  NUMBER 

8809 


5e.  TASK  NUMBER 

PPM000 18802 


5f.  WORK  UNIT  NUMBER 

EF 122707 


8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 

AFRL/RVSV 


11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 

AFRL-RV-PS-TR-20 15-0182 


14.  ABSTRACT 

Characterize  the  nonlinear  dynamics  for  large  amplitude  relative  motion  in  the  Hill  frame  (or  curvilinear  frame)  in  terms  of 
Nonlinear  Normal  Modes  for  the  case  of  a  circular  chief  orbit.  Characterize  the  nonlinear  dynamics  for  nonlinear 
translationalrotational  coupling  of  relative  motion  in  terms  of  Nonlinear  Normal  Modes.  Characterize  the  nonlinear  dynamics 
for  an  elliptic  Keplerian  chief  orbit  in  terms  of  time  periodic  Nonlinear  Normal  Modes  using  the  Liapunov-Floquet 
transformation.  If  time  permits,  investigate  the  relative  navigation  problem  in  the  Nonlinear  Normal  Modes  NNMs 
framework,  using  the  strategies  outlined  above. 


15.  SUBJECT  TERMS 

Nonlinear  Normal  Modes;  NNMs;  Tschauner-Hempel-Lawden;  THL;  Relative  Orbit  Elements;  ROEs;  Linear  Normal 
Modes;  Degree-of-freedom  Vibratory  Systems;  Liapunov-Floquet  transformation;  LFT 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT 

18.  NUMBER 
OF  PAGES 

19a.  NAME  OF  RESPONSIBLE  PERSON 

Thomas  A.  Lovell 

a.  REPORT 

Unclassified 

b.  ABSTRACT 

Unclassified 

c.  THIS  PAGE 

Unclassified 

Unlimited 

50 

19b.  TELEPHONE  NUMBER  (include  area 
code) 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  239.18 


This  page  is  intentionally  left  blank 


Approved  for  public  release;  distribution  is  unlimited. 


TABLE  OF  CONTENTS 


Section  Page 

List  of  Figures . ii 

List  of  Tables . iii 

1.0  SUMMARY . 1 

2.0  INTRODUCTION . 2 

3.0  METHODS,  ASSUMPTIONS,  AND  PROCEDURES . 3 

3.1  Dynamics  of  Relative  Motion  with  Different  Nonlinearities . 3 

3.2  Measurements . 4 

3.3  Analytical  Observability  Criteria  for  Nonlinear  Systems . 5 

3.4  Numerical  Observability  Measure  of  Nonlinear  Systems . 6 

4.0  OBSERVABILITY  ANALYSIS  FOR  RELATIVE  MOTION  WITH  LOS 
MEASUREMENTS . 8 

5 .0  NUMERICAL  RESULTS  FOR  OBSERVABILITY  ANALYSIS  OF  RELATIVE 

MOTION . 10 

5.1  Results  for  Four  Models  with  Different  Nonlinearities . 11 

5.2  Effects  of  Nonlinearities  on  the  Observability  of  Different  Configurations . 14 

6.0  AMBIGUOUS  ORBITS  OF  HCW  MODEL  WITH  RANGE-ONLY 
MEASUREMENTS . 18 

6.1  Definition  of  Ambiguous  Orbits . 18 

6.2  Mirror  Ambiguous  Orbits . 21 

6.3  Deformed  Ambiguous  Orbits . 23 

6.4  Existence  of  Deformed  Ambiguous  Orbits . 27 

6.5  Categorization  of  Ambiguous  Relative  Orbits . 29 

6.6  Numerical  Results  for  Ambiguity  Analysis . 31 

7.0  RESULTS  AND  DISCUSSION . 35 

8.0  CONCLUSIONS . 38 

REFERENCES . 39 

LIST  OF  SYMBOLS,  ABBREVIATIONS,  AND  ACRONYMS . 40 


Approved  for  public  release;  distribution  is  unlimited. 


LIST  OF  FIGURES 


Figure  Page 

Figure  1.  Measurement  Relationship  in  Chief’s  LVLH  Frame . 5 

Figure  2.  The  True  Relative  Orbit  in  Chief’s  LVLH  Frame  with  3-D  Projections . 10 

Figure  3.  True  and  Estimated  Orbits  with  Range-only  Measurements . 12 

Figure  4.  True  and  Estimated  Orbits  with  Angle s-only  Measurements . 13 

Figure  5.  Ratio  Plots  of  01  and  CN  with  Varying  Inclination  Difference . 14 

Figure  6.  Ratio  Plots  of  01  and  CN  with  Varying  Mean  Anomaly  Difference . 15 

Figure  7.  Ratio  Plots  of  01  and  CN  with  Varying  Eccentricity  Difference . 16 

Figure  8.  Plots  of  01  and  CN  for  Different  Orders  of  Nonlinearities  (Range-only) . 16 

Figure  9.  Plots  of  01  and  CN  for  Different  Orders  of  Nonlinearities  (Angles-only) . 17 

Figure  10.  Maximum  Relative  Distance  in  Three  Directions  with  Increasing  Si . 17 

Figure  11.  First  Case  of  Ambiguity  Using  HCW  Based  EKF . 18 

Figure  12.  Second  Case  of  Ambiguity  Using  HCW  Based  EKF . 19 

Figure  13.  Third  Case  of  Ambiguity  Using  HCW  Based  EKF . 20 

Figure  14.  True  Relative  Orbit  and  Three  Mirror  Ambiguous  Orbits . 23 

Figure  15.  True  Relative  Orbit  and  Four  Deformed  Ambiguous  Orbits . 26 

Figure  16.  Three  Polynomials  and  Corresponding  Solutions . 27 

Figure  17.  Cases  of  Transverse  and  Tangent  Intersections  for  h(K)  and  c . 28 

Figure  18.  Illustration  of  Tangent  Condition . 29 

Figure  19.  Illustration  of  Offset  yd0  and  Slope  s  in  y  —  z  Projection . 30 

Figure  20.  Drifting  True  and  Mirror  Ambiguous  Orbits  in  an  EKF  Simulation . 31 

Figure  21.  Type  (c)  Mirror  Ambiguous  Orbit  in  an  EKF  Simulation . 32 

Figure  22.  Type  (e)  Deformed  Ambiguous  Orbit  in  an  EKF  Simulation . 33 

Figure  23.  Relative  Orbits  Resulting  from  Variation  of  Initial  Conditions . 34 

Figure  24.  Convergence  of  EKF  Simulations  for  Small  Relative  Orbit  Scenario . 35 

Figure  25.  Convergence  of  EKF  for  Medium  and  Farge  Relative  Orbit  Scenarios . 36 

Figure  26.  Convergence  of  EKF  for  HCW  Model  with  Three  Relative  Orbit  Scenarios . 36 

Figure  27.  Convergence  of  EKF  for  Full  Nonlinear  Model  with  Three  Scenarios . 37 


Approved  for  public  release;  distribution  is  unlimited. 

ii 


LIST  OF  TABLES 


Table  Page 

Table  1.  Orbital  Element  Differences  for  Primary  Simulation  Case . 10 

Table  2.  Performance  of  EKF  Based  on  Range-only  Measurements . 11 

Table  3.  Performance  of  EKF  Based  on  Angles-only  Measurements . 13 

Table  4.  Types  of  Mirror  Ambiguous  Orbits  in  Initial  Cartesian  Coordinates . 22 

Table  5.  Initial  Condition  for  True  Relative  Orbit . 23 

Table  6.  Categorization  of  Ambiguous  Orbits . 30 

Table  7.  Orbital  Element  Differences  for  Drifting  and  Non-drifting  Scenarios . 32 

Table  8.  Orbital  Element  Differences  for  Three  Simulation  Scenarios . 34 


Approved  for  public  release;  distribution  is  unlimited. 

iii 


ACKNOWLEDGMENTS 


This  material  is  based  on  research  sponsored  by  Air  Force  Research  Laboratory  under  agreement 
number  FA9453-15-1-0350.  The  U.S.  Government  is  authorized  to  reproduce  and  distribute  reprints 
for  Governmental  purposes  notwithstanding  any  copyright  notation  thereon. 

DISCLAIMER 

The  views  and  conclusions  contained  herein  are  those  of  the  authors  and  should  not  be  interpreted  as 
necessarily  representing  the  official  policies  or  endorsements,  either  expressed  or  implied,  of  Air 
Force  Research  Laboratory  or  the  U.S.  Government. 


Approved  for  public  release;  distribution  is  unlimited. 


1.0  SUMMARY 


In  this  report,  the  effects  of  incorporating  nonlinearities  in  sequential  relative  orbit  esti¬ 
mation  are  studied  for  a  chief  spacecraft  in  a  circular  orbit,  assuming  either  range-only  or 
angles-only  measurement  of  the  deputy  from  the  chief.  The  relative  motion  models  can  be 
categorized  into  four  cases:  first  order  Hill-Clohessy-Wiltshire  (HCW)  equation,  second 
order,  third  order  and  full  nonlinear.  Observability  is  studied  analytically  using  Lie  deriva¬ 
tives  and  numerically  with  the  observability  index  and  condition  number  obtained  from 
employing  an  Extended  Kalman  Filter  (EKF).  The  results  highlight  the  improving  benefits 
of  using  higher  order  nonlinear  models. 

To  explain  the  behavior  of  HCW  dynamics  in  an  EKF  with  range-only  measurements,  a 
classification  of  ambiguous  spacecraft  relative  orbits  in  sequential  orbit  estimation  is  for¬ 
mulated  based  on  continuous  range-only  measurements.  Using  relative  orbit  elements  the 
ambiguous  orbits  are  categorized  into  two  cases:  mirror  ambiguous  orbits,  which  conserve 
the  size  and  shape  but  transform  the  orientation  of  the  true  relative  orbit,  and  deformed 
ambiguous  orbits,  which  both  distort  the  shape  and  change  the  orientation.  Furthermore,  it 
is  shown  that  the  inclusion  of  higher  order  nonlinearities  in  the  filter  model  can  help  resist 
the  tendency  of  an  EKF  to  converge  to  the  ambiguous  orbits. 
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2.0  INTRODUCTION 


Relative  orbit  estimation  is  desirable  for  many  types  of  spacecraft  missions,  such  as  for¬ 
mation  control  and  rendezvous.  Performing  spacecraft  maneuvers  based  only  on  on-board 
measurements  reduces  the  total  operating  cost,  and  improves  safety  against  communica¬ 
tion  interruptions  with  ground  stations.  Relative  navigation  between  spacecraft  in  close- 
proximity  essentially  corresponds  to  space-based  orbit  determination.  In  particular,  range- 
based  and  vision-based  navigation  and  estimation  of  relative  orbit  have  received  attention 
recently,  since  they  have  desirable  properties  of  low  cost  and  minimal  maintenance. 

Several  papers  have  considered  relative  orbit  navigation  and  estimation  from  different 
perspectives.  The  core  problem  is  to  determine  the  relative  orbit  between  a  chief  spacecraft 
and  a  deputy  spacecraft  by  certain  measurements,  assuming  that  the  orbit  of  the  chief  is 
prescribed  exactly.  Huxel  and  Bishop  [1]  discussed  the  effects  of  using  both  inertial  range 
measurements  from  tracking  stations  and  relative  range  measurements  among  formation 
members  in  the  context  of  two-body  dynamics  in  the  inertial  frame.  Also  using  two-body 
inertial  dynamics,  Yim,  et  al.  [2]  numerically  studied  the  observability  of  relative 
orbit  estimation  by  taking  line-of-sight  (LOS)  measurements  with  incorporation  of  J2 
perturba-tion  and  showed  that  proper  choices  of  orbital  element  differences  can  improve 
estimation  performance.  Taking  LOS  measurement  only,  Woffiden  and  Geller  [3,4] 
discussed  relative  orbit  estimation  based  on  the  Hill-Clohessey-Wiltshire  dynamic  model 
[5,6]  and  concluded  this  scenario  is  unobservable.  Using  the  Lie  derivative  method, 
Kaufman,  et  al.  [7]  showed  that  with  LOS  measurement  only,  the  nonlinear  relative 
orbital  dynamics  are  observable  under  certain  geometric  conditions.  Psiaki  [8]  considered 
relative  orbit  estimation  from  the  view  of  orbital  element  differences  and  widely  discussed 
the  improvement  of  observability  by  adding  J2  perturbation  or  altering  orbital  element 
differences.  Rundberg  and  Lovell  [9]  discussed  the  initial  relative  orbit  determination  using 
minimal  number  of  range-only  measurements. 

In  this  report,  the  effects  of  including  nonlinearities  in  the  filter  dynamic  model  on  ob¬ 
servability  in  relative  orbit  estimation  for  unperturbed  circular  chief  orbits  are  explored. 
Four  different  dynamic  models,  i.e.,  first  order  (HCW),  second  order,  third  order  and  full 
nonlinear  models  are  employed  in  an  extended  Kalman  filter  along  with  two  different  types 
of  measurements:  range-only  and  angles-only.  The  analytical  method  of  Lie  derivatives  and 
the  numerical  methods  of  observability  index  and  condition  number  are  applied  to  analyze 
the  observability  in  relative  orbit  estimation  with  the  four  different  models  listed  above. 

To  explain  the  appearance  of  ambiguous  trajectories  using  continuous  range-only  mea¬ 
surements  in  an  EKF,  an  analytical  analysis  of  ambiguous  conditions  is  presented  using 
relative  orbital  elements.  Subsequently,  the  enumeration  and  classification  of  these  tra¬ 
jectories  is  provided  both  by  using  Cartesian  coordinates  and  geometric  properties  of  the 
relative  orbit.  The  condition  of  existence  for  the  deformed  ambiguous  orbit  is  also  shown 
through  the  solution  of  a  fourth  order  polynomial.  Finally,  as  means  to  exclude  ambigu¬ 
ities,  we  explore  the  possibility  of  using  higher  order  nonlinear  models  to  guarantee  the 
uniqueness  of  the  estimated  orbit. 
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3.0  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 


Consider  two  satellites  orbiting  the  Earth  where  each  satellite  is  modeled  as  a  point 
mass.  Suppose  that  a  chief  spacecraft  is  on  a  circular  orbit  with  a  pre -determined  orbital 
radius  rc  G  M.  In  fact,  there  exist  several  ways  to  model  spacecraft  relative  motion,  for 
example  describing  the  motion  of  each  spacecraft  individually  in  the  inertial  frame,  or 
using  a  dynamic  model  of  orbital  element  differences.  However,  to  better  visualize  the 
relative  motion  from  the  view  of  the  chief  spacecraft,  the  relative  motion  is  described  in  the 
chief  Local  Vertical  Local  Horizontal  (LVLH)  frame  in  this  paper. 

A  LVLH  frame  is  defined  as  follows.  Its  origin  is  located  at  the  chief  satellite.  The  axaxis 
is  along  the  radial  direction  from  the  Earth  to  the  chief,  and  the  y- axis  is  along  the  velocity 
vector  of  the  chief.  The  z-axis  is  normal  to  the  orbital  plane,  and  it  is  parallel  to  the  angular 
momentum  vector  of  the  chief.  The  angular  velocity  of  the  LVLH  frame  with  respect  to 
an  inertial  frame,  expressed  in  LVLH  coordinates,  is  given  by  u>  =  [0,  0,  n]T  G  M3 ,  where 
n  =  \J /i / r 3  is  the  mean  motion  of  the  chief  satellite,  and  //  denotes  the  gravitational 
parameter  of  the  Earth.  Note  the  inertial  velocity  of  the  chief  expressed  in  the  LVLH  frame 
is  given  by  vc  =  [0,  nrc,  0]r  G  M3.  Let  the  relative  position  of  a  deputy  spacecraft  with 
respect  to  the  chief  spacecraft  be  given  by  p  =  [x,y,z]T  G  M3  in  the  LVLH  frame. 


3.1  Dynamics  of  Relative  Motion  with  Different  Nonlinearities 


In  order  to  describe  the  dynamics  of  relative  motion  concisely,  the  state  space  form 
X  =  F(X)  is  used,  where  X  =  [pT  pT]T  =  [x,y,  z,  x,y,  z]T  G  M6  are  the  position  and 
velocity  states  of  the  deputy  relative  to  the  chief  in  chief’s  LVLH  frame,  and  F(X)  is  the 
corresponding  vector  field.  The  derivation  for  dynamic  models  of  different  nonlinear  orders 
are  based  on  truncating  a  Taylor  series  expansion  of  two-body  relative  orbit  dynamics  at 
different  orders. 


First  order  dynamic  model.  The  first  order  dynamic  model  (HCW  model)  only  includes 
the  linear  terms  from  the  Taylor  expansion  in  the  relative  orbit  dynamics,  and  its  vector  field 
is  given  by 


Fi(X) 


p 

p 

_  p  _ 

—2ll>  x  p  +  K  p 

(1) 


where  K  =  diag  [  3 n2,  0,  —  n2  ]  is  a  diagonal  matrix. 

Second  order  dynamic  model.  The  second  order  dynamic  model  [10]  adds  second  order 
terms  from  the  Taylor  expansion  to  the  first  order  model,  and  its  vector  field  is  given  by 


F2(X) 


P 

P 

_  P  _ 

— 2 U)  x  p  +  K p  +  r2 

(2) 


where  T2  represents  2nd  order  terms,  namely 


r2 


if 

(v»4 


-3a;2  +  |  y1  +  \z 
3  xy 
3  xz 


2  1 


(3) 
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Third  order  dynamic  model.  The  third  order  dynamic  model  [11]  adds  third  order  terms 
from  the  Taylor  expansion  to  the  second  order  dynamic  model,  and  its  vector  field  is  given 
by 


F3(X) 


p 

p 

_  p  _ 

— 2cu  x  p  +  K p  +  T2  +  r3 

(4) 


where  T3  represents  third  order  terms 


r3 


4x3  —  6  x(y2  +  z 2) 
—6  x2y  +  +  \yz 2 

—6  x2z  +  fz3  +  \zy2 


(5) 


Full  nonlinear  dynamic  model.  The  full  nonlinear  dynamic  model  is  the  most  accurate 
model  considered  in  this  paper.  The  vector  field  can  be  expressed  as 


F/(X) 


p 

p 

_  p  _ 

[  2c a  x  p  UJ  x  (u>  x  rd)  f|^jj3  J 

(6) 


where  r,/  =  [rc  +  x,  y,  z]T  G  M3  is  the  position  vector  of  the  deputy  relative  to  the  center  of 
the  Earth  in  the  chief’s  LVLH  frame. 


3.2  Measurements 

For  the  output  equation,  we  assume  the  chief  spacecraft  takes  relative  measurements 
towards  the  deputy  of  either  range-only  measurements  or  angles-only  measurements. 

Range  measurements.  Range  measurements  are  represented  by  the  magnitude  of  the 
relative  position  vector,  i.e., 


Y  =  p  =  y/x2  +  y2  +  z2 


(7) 


Angle  measurements.  One  option  to  define  angle  measurements  is  by  using  two  bearing 
angles  A  and  (  shown  in  Figure  1),  namely  the  measurement 


Y  = 


(8) 


where 


</> 

A 

A 


z 

asm  = 

V  a:2  +  y2  + 

atan2  (y,x),  if 

27t  +  atan2 (y,  x), 


atan2 (y,x)  G  [0,7r] 
if  atan2(y,  x)  G  (— n,  0) 


(9) 
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Figure  1.  Measurement  Relationship  in  Chief’s  LVLH  Frame 

Another  option  is  the  equivalent  LOS  measurement,  which  is  the  unit- vector  in  the  rela¬ 
tive  position  direction,  i.e., 


Y  =  p  =  (10) 

IIpII 

It  is  clear  to  see  that  the  LOS  vector  p  and  two  bearing  angles  (A,  6)  are  equivalent. 

3.3  Analytical  Observability  Criteria  for  Nonlinear  Systems 

Generally,  for  a  nonlinear  system  given  by 


X 

=  f{x) 

y 

=  h(x) 

x(0) 

=  x0 

(11) 

where  x  e  M"  and  y  e  W,  the  system  is  observable  over  the  interval  [0,  T]  if  the  mapping 
from  initial  state  xq  to  output  profile  y(  0  :  T)  is  one  to  one.  It  is  locally  observable  over  the 
interval  [0,  T]  if  this  mapping  is  locally  one  to  one.  A  widely  accepted  tool  for  checking 
local  observability  is  the  observability  rank  condition  given  by  Lie  derivatives  [12, 13].  Let 
the  first  order  Lie  derivative  of  output  h(x )  along  vector  field  f(x)  be 

LMx)  =  ^^f(x)eRpxl  (12) 

with  the  nth  order  Lie  derivative  defined  recursively  as 

dUi-'Hx) 

Lnfh(x )  =  — \  fix )  G  Mpxl  (13) 
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with  zero  order  Lie-derivative  L°jh  =  h.  Define  an  observability  matrix  N  e  Wnxn 


as 


N(x) 


d_ 

dx 


L°fh(x) 

Lfh(x) 

Lnflh(x) 


(14) 


It  has  been  shown  that  the  system  is  locally  observable  at  x  if  the  rank  of  the  observabil¬ 
ity  matrix  satisfies  rank  N(x)  =  n.  When  applied  to  linear  dynamics,  this  yields  the 
well-known  observability  rank  condition  for  linear  systems.  Note  that  when  there  is  more 
than  a  single  measurement  type,  i.e.,  p  >  1,  it  is  possible  to  satisfy  the  observability  rank 
condition  without  need  for  computing  the  higher-order  Lie  derivatives  up  to  the  (n  —  l)th 
order.  Moreover,  it  is  usually  difficult  to  check  higher  order  rows  in  N (x),  and  the  practical 
choice  is  to  stop  when  N(x o)  has  more  rows  than  columns,  or  at  least  when  N(x)  is  a 
square  matrix.  For  this  reason,  the  Lie  derivative  method  can  only  guarantee  observability 
if  the  first  n  rows  in  matrix  N (x)  are  non-singular,  but  the  method  cannot  guarantee  the 
unobservability  of  the  nonlinear  system. 


3.4  Numerical  Observability  Measure  of  Nonlinear  Systems 

The  observability  rank  condition  essentially  determines  whether  the  system  is  locally 
observable  over  a  short  period  of  time  in  the  vicinity  of  the  time  at  which  the  observability 
rank  is  computed,  but  it  does  not  tell  us  how  easy  it  is  to  observe  the  system.  To  overcome 
this  problem,  the  observability  Gramian  and  some  related  quantities  are  introduced. 

The  observability  Gramian  measures  the  sensitivity  of  the  output  with  respect  to  the  ini¬ 
tial  condition.  For  a  continuous  nonlinear  system  as  in  Eq.  (11),  the  observability  Gramian 
is  defined  as  [14]: 


dy(r) 

dx0 


dr 


(15) 


Corresponding  to  this  definition,  the  observability  Gramian  for  discrete  time  systems  is 
defined  as: 


N 

Wd(X0,ti,tf)  =  (16) 

i= 0 

where  N  is  the  number  of  measurement  times,  <F(i,  0)  is  the  state  transition  matrix  from 
time  t  —  to  to  time  t  —  tt  satisfying  6(f)  =  ^^  <&(£),  and  H  —  is  the  Jacobian  matrix 
for  measurement  relationship.  It  turns  out  that  the  observability  of  a  system  is  deeply 
related  to  the  eigenvalues  and  singular  values  of  the  observability  Gramian,  with  which 
a  wide  variety  of  different  measures  have  been  proposed.  Two  of  most  commonly  used 
measures  are  introduced  here. 
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Waldraff,  et  al.  [15]  outlined  the  observability  index  (01),  which  is  defined  as  the  small¬ 
est  singular  value  of  the  observability  Gramian  Wd,  i.e., 

01  =  min  a(Wd)  (17) 

where  a  denotes  singular  value  of  a  matrix.  They  also  discussed  that  if  01  is  small,  then 
observation  noises  can  have  a  large  impact  on  the  estimation  error.  In  other  words,  a  larger 
01  indicates  better  observability. 

Dochain,  et  al.  [16]  made  use  of  the  estimation  condition  number  (CN)  for  the  ob¬ 
servability  analysis,  which  is  defined  as  the  ratio  of  largest  singular  value  to  the  smallest 
singular  value  of  the  observability  Gramian  Wd,  i.e., 

CN= ZZfj™  (18) 

nun  a(Wd) 

If  CN  is  large  then  the  effect  on  the  output  caused  by  a  small  change  in  the  initial  condition 
in  one  direction  can  swamp  the  effect  on  the  output  of  a  change  in  another  direction.  In 
other  words  the  estimation  problem  is  ill-conditioned  near  states  with  large  local  estimation 
condition  number. 
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4.0  OBSERVABILITY  ANALYSIS  LOR  RELATIVE  MOTION  WITH  LOS  MEA¬ 
SUREMENTS 


As  discussed  in  the  Introduction,  Reference  [7]  applied  the  Lie  derivative  method  on  the 
full  nonlinear  dynamic  model,  resulting  in  observability  conditions  for  LOS  measurements. 
This  paper  extends  this  work  to  three  other  dynamic  models  and  discusses  observability 
conditions  for  each  one  of  them.  The  general  logic  is  shown  as  follows: 

Consider  a  general  form  of  vector  field  F 


F(X) 


p 

p 

_  p  _ 

— 2cu  x  p  +  . . . 

(19) 


For  the  relative  motion  models  with  different  orders  of  nonlinearities,  F  can  be  replaced  by 
Fi,  F2,  F3  or  F  f,  and  the  LOS  measurement  Y  can  be  expressed  as 

Y  =  p  =  TTTi  (20) 

IIpII 


With  certain  order  of  vector  field  F  and  LOS  measurement  Y,  we  can  derive  the  first  three 
rows  of  N  by 


-  dp 

dp  ' 

dp 

dp 

dp 

dp 

% 

dp 

dp 

dp 

-  dp 

dp  . 

The  key  to  observability  analysis  is  to  find  the  conditions  that  guarantee  the  full  rank  of  the 
N3  matrix.  However,  the  process  of  deriving  the  full  rank  condition  for  N:i  is  complicated 
and  therefore  not  shown  here.  Readers  may  refer  to  Reference  [7]  to  have  a  better  under¬ 
standing.  Instead,  we  extend  the  observability  conditions  to  relative  motion  models  with 
different  orders  of  nonlinearities. 

Obsen’ability  conditions  for  first  order  model.  The  observability  problem  with  linear 
dynamics  based  on  LOS  measurements  has  been  widely  discussed.  Most  notably,  Woffiden 
and  Geller  [3]  have  proved  that  based  on  linear  dynamics  of  the  HCW  model,  the  system 
is  unobservable  by  taking  angles-only  measurements.  In  this  paper,  by  substituting  vector 
field  Fi  into  Eq.  (21),  we  find  that  the  N:i  matrix  loses  full  rank  for  linear  dynamics, 
meaning  that  the  Lie  derivative  analysis  also  implies  the  unobservability  of  linear  dynamics 
with  angles-only  measurements. 

Obsen’ability  conditions  for  second  order  model.  Substituting  vector  field  F2  in  Eq.  (3) 
into  Eq.  (21),  the  resulting  observable  conditions  are  obtained  as: 

(i)  when  p  x  p  =  0,  Y-2^  p  and  pT{Y2  x  (u)  x  p))  0 

(ii)  when  p  x  p  f  0,  T2)^  p  and  pT( r2  x  vrei)  f  0  (22) 

where  vre/  =  p  +  u>  x  p  is  the  velocity  vector  of  deputy  relative  to  chief. 
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Obsen’ability  conditions  for  third  order  model.  Substituting  vector  field  F3  in  Eq.  (4) 
into  Eq.  (21),  the  resulting  observable  conditions  are  obtained  as: 

(i)  when  p  x  p  =  0,  (T2  +  2r3)  ft  p  and  pT(( V2  +  2r3)  x  (a;  x  p))  f  0 

(ii)  when  p  x  p  f  0,  (T2  +  2r3)  ft  p  and  pT(( V2  +  2r3)  x  vret)  f  0  (23) 


Obsen’ability  conditions  for  full  nonlinear  model.  Substituting  vector  field  F  f  in  Eq. 
(6)  into  Eq.  (21),  the  resulting  observable  conditions  are  obtained  as: 

( i )  when  p  x  p  =  0,  aj  ft  p  and  pT{cHf  x  (w  x  p))  f  0 

(ii)  when  p  x  p  f  0,  a/  jj'  p  and  pT(a/  x  vre;)  f  0  (24) 


where  a/  can  be  expressed  as 


af  =  lo  x 


(cu  x  rc)  +  ^  + 


prc  3prdrqa 


(25) 


It  is  noted  that  even  if  Eq.  (22),  (23)  or  (24)  is  violated,  it  only  implies  that  the  correspond¬ 
ing  order  of  dynamic  model  is  unobservable  at  that  specific  time  epoch  and  does  not  mean 
the  whole  measurement  profile  can  not  result  in  an  observable  system.  In  other  words, 
when  the  observable  conditions  are  violated  at  a  certain  time  epoch,  the  LOS  measure¬ 
ments  make  little  or  no  contribution  on  improving  the  observability.  On  the  contrary,  for 
linear  dynamics,  Woffiden  and  Geller  [3]  have  drawn  a  decisive  conclusion  that  the  system 
is  unobservable  with  LOS  measurements  regardless  of  the  number  of  measurements  being 
taken. 
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5.0  NUMERICAL  RESULTS  FOR  OBSERVABILITY  ANALYSIS  OF  RELATIVE 
MOTION 

We  choose  one  primary  scenario  of  relative  motion  between  chief  and  deputy  orbits 
shown  in  Table  1.  Note  that  the  orbital  elements  difference  are  chosen  such  that  the  scales 
of  relative  motion  on  three  directions  of  chief’s  LVLH  frame  are  comparable  as  shown  in 
Figure  2. 


Table  1.  Orbital  Element  Differences  for  Primary  Simulation  Case 


a  (km) 

e 

i 

it 

UJ 

M r 

Chief 

7100 

0 

0° 

0° 

0° 

0° 

Deputy 

7100 

0.005 

0.3° 

0° 

0° 

0.1° 

x  10 


x  10 


X  (m) 


x  10 


x(m) 


x  10 


y  (m) 


x  10 


x  10 


Figure  2.  The  True  Relative  Orbit  in  Chief’s  LVLH  Frame  with  3-D  Projections 


Furthermore,  we  evaluate  the  observability  measure  for  different  order  of  nonlinearities 
based  on  local  observability  Gramian  and  process  range-only  and  angles-only  measure¬ 
ments  with  the  EKF,  using  the  full  nonlinear  model  as  the  truth  model. 

For  both  range-only  and  angles-only  measurements,  the  initial  state  deviation  for  EKF  is 
5X0  =  [1000,  -1000, 1000,  l,l,-l](m,  m/s)  (26) 

and  the  initial  covariance  of  the  state  is  chosen  as 

P0  =  diag  [1.8  x  107, 1.8  x  107,1.8  x  107, 18, 18, 18]  (m2,  m2/s2)  (27) 
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The  sample  time  is  At  =  10  seconds.  For  range-only  measurements,  the  measurement 
covariance  matrix  is  assumed  to  be  R  =  202  m2  .  For  angles-only  measurements,  R  = 
diag  [  2.7416  x  10-',  2.7416  x  10~7  ]  rad2,  which  corresponds  to  the  variances  of  two 
bearing  angles  <J\  =  =  0.03°.  Since  different  orders  of  nonlinearities  will  be  considered, 

different  levels  of  process  noise  will  be  applied  to  different  dynamical  models  separately. 


5.1  Results  for  Four  Models  with  Different  Nonlinearities 


Figure  3  illustrates  the  simulation  results  for  the  EKF  based  on  the  four  different  rela¬ 
tive  motion  models  with  range-only  measurements,  in  which  the  red  line  denotes  the  true 
relative  orbit  and  the  blue  line  denotes  the  estimated  orbit.  From  Figure  3(a),  it  is  clear 
that  the  Kalman  filter  based  on  the  HCW  model  cannot  estimate  the  true  states  effectively 
and  manifests  the  inability  to  distinguish  between  two  orbits  with  the  same  range  informa¬ 
tion.  In  Figure  3(b),  using  the  second  order  dynamic  model,  the  filter  successfully  captures 
the  motion  of  true  relative  orbit,  even  though  its  accuracy  may  not  be  satisfactory  at  the 
beginning.  In  Figure  3(c),  with  a  dynamic  model  truncated  at  third  order,  the  estimated 
relative  orbit  very  closely  follows  the  true  relative  orbit.  In  fact,  it  is  surprising  that  the 
estimated  results  in  (c)  are  almost  as  good  as  the  results  in  (d).  This  indicates  that  a  third 
order  model  is  required  to  accurately  characterize  the  relative  motion  for  this  particular 
case.  It  is  expected  that  for  cases  where  the  chief  and  deputy  are  closer  together,  the  results 
utilizing  a  second  order  model  may  be  more  similar  to  those  obtained  with  a  third  order  or 
full  nonlinear  model. 


Table  2  summarizes  the  performance  of  the  EKF  with  different  nonlinearities  for  range- 
only  measurements.  01  and  CN  are  the  observability  index  and  condition  number  of  the 
observability  Gramian.  ep  indicates  the  estimation  position  error  root  mean  square  (RMS) 
and  is  defined  as  . _ 


ep  = 


\ 


N 


N 


£ll(p*-P<0ll: 


k= 0 


(28) 


where  pk  =  [xk,  Vk,  4 a-]7'  is  the  Ath  estimated  relative  position  vector  and  pk  is  the  A  th  true 
relative  position  vector.  ep  denotes  the  estimation  velocity  error,  i.e., 


eA  = 


\ 


N 


N 


£||(Ac-A0II; 


fc=0 


where  pk  =  [xk,  yk,  zk] "J 


(29) 


Table  2.  Performance  of  EKF  Based  on  Range-only  Measurements 


01 

CN 

ep  (m) 

ep  (m/s) 

^ max  (P) 

det(P) 

2  st 

1.89  x  10“i5 

2.93  x  10iS 

5.15  x  104 

54.20 

2.36  x  103 

6.00  x  10"4 

2nd 

1.61  x  10-12 

1.38  x  1013 

1.11  x  102 

0.62 

4.73  x  102 

6.16  x  10“9 

3rd 

1.33  x  io~9 

3.78  x  1010 

23.26 

0.064 

8.63 

2.41  x  10“19 

full 

1.63  x  10“9 

3.43  x  1010 

21.83 

0.052 

7.62 

5.23  x  IQ"20 
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In  Table  2,  A mai,(P)  is  the  maximum  eigenvalue  of  covariance  matrix  P.  It  indicates 
the  variance  of  worst  estimated  state  among  all  the  six  states.  det(P)  is  the  determinant 
of  matrix  P,  which  can  be  interpreted  as  the  overall  performance  of  EKF  since  det(P)  = 
A1A2  •  •  •  A6  is  the  product  of  the  variances  of  the  six  states. 

From  the  first  two  columns  of  Table  2,  01  increases  and  CN  decreases  drastically  from 
first  order  to  second  order  and  second  order  to  third  order  dynamic  models.  However, 
changing  from  third  order  model  to  full  nonlinear  model,  01  and  CN  only  differ  slightly. 
It  is  noted  that  all  the  other  parameters,  ep,  ep,  \rnax ( P)  and  det(P)  all  follow  the  same 
pattern  as  with  01  and  CN.  This  observation  agrees  with  the  simulation  results  in  Figure 
3. 


(c)  Third  order  (d)  Full  nonlinear 

Figure  3.  True  and  Estimated  Orbits  with  Range-only  Measurements 


With  angles-only  measurements,  the  simulation  results  are  illustrated  in  Figure  4  and 
Table  3.  Again,  the  observability  changes  greatly  from  the  first  order  to  second  order  and 
from  second  order  to  third  order  dynamic  models,  while  the  estimation  results  of  third  order 
and  full  nonlinear  order  models  are  almost  indistinguishable. 


Approved  for  public  release;  distribution  is  unlimited. 

12 


However,  compared  with  the  range-only  results  for  the  HCW  model  shown  in  Figure 
3(a),  the  estimated  orbits  with  angles-only  measurements  illustrated  in  Figure  4(a)  are 
shrunk  and  keep  the  shape  of  the  true  orbit  instead  of  being  disoriented.  This  is  mainly 
because  of  the  incapability  of  EKF  with  angles-only  measurement  on  capturing  range  in¬ 
formation,  which  has  been  proposed  by  Woffiden  and  Geller  [3]  as  a  scalar  ambiguity  of 
angles-only  measurements  based  on  the  linear  dynamics. 


(c)  Third  order  (d)  Full  nonlinear 

Figure  4.  True  and  Estimated  Orbits  with  Angles-only  Measurements 


Table  3.  Performance  of  EKF  Based  on  Angles-only  Measurements 


01 

CN 

ep  (m) 

ep  (m/s) 

^max  (P ) 

det(P) 

2  st 

2.27 

X 

O 

T— i 

=rTd 

5.24 

X 

10iv 

4.27  x  104 

44.47 

781.43 

1.83  x  10~9 

2nd 

1.75 

X 

10 

-14 

1.41 

X 

1014 

1.96  x  102 

1.95 

136.51 

2.48  x  10"13 

3rd 

7.32 

X 

10 

-10 

3.85 

X 

1010 

3.52 

0.035 

42.735 

2.27  x  10-23 

full 

7.41 

X 

10 

-10 

3.62 

X 

1010 

2.81 

0.029 

38.729 

3.40  x  10-24 
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5.2  Effects  of  Nonlinearities  on  the  Observability  of  Different  Configurations 

In  this  section,  we  study  how  observability  changes  with  different  nonlinear  filter  models 
through  the  variation  of  orbital  element  differences.  The  methodology  is  to  vary  a  certain 
orbit  element  difference  and  compare  the  improvement  of  01  and  CN  between  adjacent 
dynamic  models  (e.g.  first  and  second  order,  second  and  third  order,  and  third  order  and 
full  nonlinear). 

Figure  5  illustrates  the  change  in  observability  when  varying  inclination  difference.  Fig¬ 
ures  5(a)  and  5(b)  are  cases  for  range  measurements,  while  Figures  5(c)  and  5(d)  are  cases 
for  angle  measurements.  From  these  plots,  we  find  that  with  increased  inclination  differ¬ 
ence  5i,  the  observability  difference  improves  almost  uniformly  (indicated  by  the  increasing 
log  Olj/OIi  and  log  CNi/CNj,  j  >  i).  This  means  that  higher  order  nonlinearities  grant 
extra  benefits  when  inclination  difference  is  enlarged.  Figure  6  and  7  also  illustrate  cases 
for  variation  of  mean  anomaly  differences  and  eccentricity  differences.  From  the  results, 
the  same  conclusion  can  be  drawn  as  with  variation  of  inclination  difference. 


8  i  (degree) 

(a)  Ratio  plots  of  OI  for  range-only 


8  i  (degree) 

(b)  Ratio  plots  of  CN  for  range-only 
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(c)  Ratio  plots  of  OI  for  angles-only  (d)  Ratio  plots  of  CN  for  angles-only 

Figure  5.  Ratio  Plots  of  OI  and  CN  with  Varying  Inclination  Difference 
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(a)  Ratio  plots  of  OI  for  range-only  measurements  (b)  Ratio  plots  of  CN  for  range-only  measure¬ 
ments 


(c)  Ratio  plots  of  OI  for  angles-only  (d)  Ratio  plots  of  CN  for  angles-only  measure¬ 

ments 

Figure  6.  Ratio  Plots  of  OI  and  CN  with  Varying  Mean  Anomaly  Difference 


Figures  8  and  9  present  the  change  of  OI  and  CN  for  variation  of  inclination  difference 
with  range-only  and  angles-only  measurements  separately.  From  these  results,  the  OI  and 
CN  of  the  third  order  model  follow  very  closely  to  those  of  the  full  nonlinear  model  for 
both  range-only  and  angles-only  measurements  as  the  inclination  difference  Si  increases. 
With  increasing  inclination  difference,  the  advantages  of  using  higer  order  models  become 
more  obvious,  particularly  for  including  second  and  third  order  terms. 

However,  the  plots  of  the  second  order  model  drift  away  from  the  full  nonlinear  model, 
meaning  that  the  second  order  model  cannot  accurately  capture  all  the  properties  of  relative 
orbital  motion  of  the  true  model.  With  no  surprise,  the  first  order  dynamic  model  gives  the 
worst  performance  among  all  the  relative  motion  models  because  of  its  assumption  of  close 
proximity  on  which  the  linearization  is  valid. 

We  also  notice  that  there  exists  tuning  point  in  Figures  8  and  9.  For  example,  in  Fig¬ 
ure  8(a),  the  maximum  value  of  OI  appears  between  Si  =  0.2°  and  Si  =  0.8°,  and  the 
corresponding  value  of  Si  differs  slightly  between  systems  with  different  nonlinearities. 
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(a)  Ratio  plots  of  OI  for  range-only  measurements 


(b)  Ratio  plots  of  CN  for  range-only  measurements 


8  e  8  e 


(c)  Ratio  plots  of  OI  for  angles-only  measurements  (d)  Ratio  plots  of  CN  for  angles-only  measurements 

Figure  7.  Ratio  Plots  of  OI  and  CN  with  Varying  Eccentricity  Difference 


(a)  OI  for  range-only  measurements  (b)  CN  for  range-only  measurements 

Figure  8.  Plots  of  OI  and  CN  for  Different  Orders  of  Nonlinearities  (Range-only) 
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To  explain  this  behavior,  recall  that  01  measures  the  worst  observability  among  all  six 
states  in  relative  orbit  estimation  and  the  filter  actually  uses  the  change  of  range  information 
to  decide  the  deputy’s  orbit,  i.e.,  the  x,  y,  z  coordinates  of  deputy  in  chief’s  LVLH  frame. 
From  p  =  \Jx2  -\-  y2  +  z 2,  we  have  bp  =  LpSe,  where  e  denotes  x,  y  or  2.  Therefore,  it  is 
clear  that  when  e  is  the  smallest  among  all  three  coordinates,  it  will  be  the  most  insensitive 
state  regarding  the  change  of  range,  and  therefore  the  least  observable  state.  From  Figure 
10,  we  can  see  that  when  5i  <  0.3°,  the  maximum  scale  in  z  direction  is  less  than  that  in 
x  and  y  direction,  making  z  coordinates  least  observable.  However,  when  Si  >  0.8°,  the 
scale  of  z  coordinates  becomes  the  dominant,  making  x  coordinates  least  observable.  This 
observation  agrees  with  the  occurrence  of  the  maximum  value  of  01  between  Si  =  0.2° 
and  Si  =  0.8°.  Furthermore,  the  slightly  different  peak  values  of  Si  come  from  the  use  of 
dynamic  models  with  different  orders  of  nonlinearities. 


(a)  OI  for  angles-only  measurements  (b)  CN  for  range-only  measurements 

Figure  9.  Plots  of  OI  and  CN  for  Different  Orders  of  Nonlinearities  (Angles-only) 


Figure  10.  Maximum  Relative  Distance  in  Three  Directions  with  Increasing  Si 
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6.0  AMBIGUOUS  ORBITS  OF  HCW  MODEL  WITH  RANGE-ONLY  MEASURE¬ 
MENTS 


6.1  Definition  of  Ambiguous  Orbits 

In  the  previous  section,  we  get  a  glimpse  of  the  weakness  of  the  HCW-based  EKF  with 
range-only  measurements  on  relative  orbit  estimation  in  Figure  3(a).  However,  when  exe¬ 
cuting  a  series  of  simulations  with  HCW  dynamic  model,  the  estimated  orbits  converge  to 
several  different  orbits  randomly. 

First  Case:  True  relative  orbit.  The  results  for  the  first  case  are  shown  in  Figure  11. 
From  Figure  11(a),  we  can  see  that  the  estimated  orbit  follows  fairly  closely  to  the  true 
orbit.  From  Figure  11(b),  even  though  the  estimation  error  is  still  in  the  scale  of  500  m,  the 
filter  successfully  keeps  the  estimated  trajectory  from  diverging. 


(a)  True  vs.  estimated  orbits  (red:  true  orbit,  blue:  (b)  Position  error  (blue:  estimation  error,  red:  3cr 
estimated  orbit)  covariance  envelop) 

Figure  11.  First  Case  of  Ambiguity  Using  HCW  Based  EKF 


Second  Case:  Mirror  ambiguous  orbit.  From  Figure  12(a),  the  estimated  orbit  (blue 
line)  looks  like  a  mirror  image  of  the  true  orbit  (red  line),  which  indicates  the  two  orbits 
have  same  shape  and  size.  Furthermore,  in  Figure  12  (b),  (c)  and  (d),  the  in-plane  xy  motion 
is  described  well  by  the  estimated  orbit,  but  the  out-of-plane  z  motion  is  poorly  estimated. 

Third  Case:  Deformed  ambiguous  orbit.  Figure  13  illustrates  the  trajectory  of  the  true 
orbit  (red  line)  and  the  deformed  estimated  orbit  (blue  line).  Compared  with  the  true  orbit, 
the  in-plane  xy  motion  of  the  estimated  orbit  is  reduced,  and  the  out-of-plane  z  motion 
is  enlarged.  It  appears  that  the  EKF  with  HCW  dynamics  overestimated  the  z  motion  to 
compensate  for  the  underestimated  xy  motion  to  meet  the  range  requirements.  Also,  the 
orientation  of  the  estimated  orbit  is  different  from  the  true  orbit. 

For  the  following  discussion,  we  will  discuss  the  existence  of  these  ambiguous  orbits 
and  their  effects  on  the  performance  of  a  filter.  For  convenience,  we  make  the  following 
definition.  An  orbit  is  defined  as  an  ambiguous  orbit  of  the  true  orbit  by  range  if  it  shares  the 
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x  104 


(a)  True  orbit  vs.  estimated  orbit  (b)  xy  projection 


x  104 


x  104 


(d)  yz  proejction 


Figure  12.  Second  Case  of  Ambiguity  Using  HCW  Based  EKF 


same  range  history  with  the  true  orbit.  Namely,  p'(t)  =  pit)  for  all  t  >  0,  where  pit)  and 
p'{t )  are  the  ranges  of  the  true  and  ambiguous  orbits  at  time  t  respectively.  Furthermore,  if 
the  ambiguous  orbit  conserves  the  size  and  shape  of  the  true  orbit  then  it  is  classified  as  a 
mirror  ambiguous  orbit;  otherwise,  it  is  classified  as  a  deformed  ambiguous  orbit. 

To  analyze  the  ambiguous  relative  orbits,  we  utilize  the  solution  of  the  HCW  equation 
(Eq.  (2)  expressed  in  terms  of  relative  orbit  elements  (ROEs)  [17], 

x(t)  =  -y  cos(^)  +Xd 

y{t)  =  ae  sin(/3)  +  yd 

z(t )  =  zm  cos(V0  (30) 

where  ae,  xd,  zm  are  constant  and  yd(t )  =  yd0  —  | nxdt,  (3(t)  =  /30  +  nt  and  ip(t)  =  ip0  +  nt 
are  time  dependent.  It  is  clear  that  (ae,  zm.  xd,  yd0,  l3(h  Wo)  are  six  constants  that  can  be  used 
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x  104 


(a)  True  orbit  vs.  estimated  orbit 


(b)  xy  projection 


xIO4  xIO4 


(c)  xz  projection  (d)  yz  proejction 

Figure  13.  Third  Case  of  Ambiguity  Using  HCW  Based  EKF 


to  represent  relative  orbits.  The  square  of  range  p  at  time  t  can  be  expressed  as 


p\t) 


=  x\t)+y2(t)  +  z2(t) 

5  o  1  9  9 

=  gae  +  +  yd0 


+ 


3  1 

-g al  cos(2/30)  +  -4cos(2^0) 


3  1 

-a2  sin(2/?0)  -  -z2n  sin(2^0) 


cos(2  nt) 
sin(2  nt) 


+2aeyd0  sin(/30)  cos  (nt)  +  2  aeyd0  cos(/30 )  sin  (nt) 

—3  aexd  sm(/3o)nt  cos  (nt)  —  3aexd  cos(/3o)nt  sin  (nt) 
9 

-Sxdyd07it  +  -x2dt2 


(31) 


Since  the  nine  basis  functions  in  Eq.  (31)  are  linearly  independent,  the  following  nine 
identities  must  be  satisfied  for  p'2{t )  =  p2{t)  . 
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5  d2e  +  4z2n  +  8yd0 

=  ho!2  +  4z'm2  +  8y'd02 

(32) 

3 a2  cos(2/30)  -  4 z2m  cos(2^0) 

=  3 a'2  cos(2/3q)  -  4 z'J  cos(2^) 

(33) 

3 a2e  sm(2f3o)  -  4 z2n  sin(2-^0) 

=  3a'e2  sin(2/3o)  -  4 z'J  sin(2^) 

(34) 

aeyd  o  sin(^o) 

=  aWdo  sin(/?o) 

(35) 

aeVd  o  cos(/30) 

=  a!ey'd  ocos(/5o) 

(36) 

aexdsm(/30) 

=  aWds  in(^o) 

(37) 

aexd  cos(/30) 

=  a'ex'd  cos(Pq) 

(38) 

•^dVdO 

=  xdyd o 

(39) 

x2d 

/  2 

=  Xd 

(40) 

where  the  primed  quantities  correspond  to  the  ambiguous  orbit.  It  is  noted  that  if  the  non¬ 
drifting  condition  xd  =  0  is  satisfied,  Eqs.  (37-40)  will  vanish.  In  the  following  section,  we 
first  discuss  the  ambiguity  under  the  non-drifting  condition  (Eqs.  (32-36))  and  then  check 
the  validity  of  the  ambiguity  once  the  non-drifting  assumption  is  violated.  First,  however, 
we  observe  that  Eqs.  (35)  and  (36)  yield 

I  aeUdO  |  =  Ky'dol  (41) 


6.2  Mirror  Ambiguous  Orbits 

A  mirror  ambiguous  orbit  must  have  the  same  magnitude  of  in-plane  and  out-of-plane 
motion  as  the  true  orbit,  i.e.,  \ae\  =  \a!e\  and  zm  =  \z'm\.  Two  cases  are  considered. 

Case  1:  aeyd0  =  a'ey'd0. 

This  case  yields  either 

a'e  =  ae 
y'do  =  Vdo 

or, 

a!e  =  —  ae 

y'do  =  -ydo 

Substituting  aeyd 0  =  a'eyd0  into  Eqs.  (35)  and  (36),  we  obtain 

Po  =  Po  ,  for  fa  e  [0, 2vr)  (46) 

Then  substituting  Eqs.  (42-43)  or  (44-45)  into  Eq.  (32) ,  we  obtain  z2m  =  z'm2  yielding 

Z'm  =  Zm  (47) 

z'm  =  -zm  (48) 


(44) 

(45) 


(42) 

(43) 
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Finally,  with  the  substitution  Eqs.  (42-48)  into  Eq.  (33),  we  have 

V’o  =  A  (49) 

V’o  =  V’o  +  7T  (50) 

for  -0  e  [0,  2n).  This  yields  eight  non-unique  possible  combinations  of  the  initial  ROEs. 


Case  2:  aeyd0  =  -a'ey'd0. 

As  in  Case  1,  this  also  yields  eight  non-unique  initial  ROE  combinations. 

These  16  cases  contain  redundancy  in  terms  of  initial  Cartesian  coordinates,  because 
the  initial  Cartesian  coordinates  uniquely  determine  an  orbit,  whereas  ROEs  do  not.  The 
transformation  from  initial  ROEs  to  initial  Cartesian  states  [17]  is 


x0  = 

-y  COS  A) 

Vo  = 

ae  sin  f30  +  yd0 

z0  = 

Zm  COS 

x0  = 

ysm  A) 

yo  = 

nae  cos  f3o 

i0  = 

—nzm  sin 

With  the  transformation  in  Eq  (51),  it  can  be  shown  that  the  total  number  of  Cartesian 
initial  states  that  result  in  orbits  with  the  same  range  history  is  four  including  the  true  orbit 
and  three  mirror  ambiguous  orbits.  Table  4  lists  these  orbits  as  (a),  ( b ),  (c)  and  (t)  (which 
stands  for  the  true  orbit).  For  the  case  of  drifting  orbit,  it  is  easy  to  check  that  the  foregoing 
16  ROEs  combinations  in  Case  1  and  2  also  satisfy  Eqs.  (37-40).  Therefore,  the  drifting 
phenomenon  does  not  exclude  these  three  mirror  ambiguous  orbits. 


Table  4.  Types  of  Mirror  Ambiguous  Orbits  in  Initial  Cartesian  Coordinates 


xo 

Xq 

y'o 

y'o 

z'o 

zo 

(0 

Xo 

Xo 

yo 

yo 

Zo 

Z  0 

(a) 

x0 

x0 

yo 

yo 

-Zq 

-Zo 

(b) 

-x0 

-x0 

-yo 

-i/o 

Zo 

Zo 

(c) 

-x0 

-x0 

-yo 

-yo 

-Zo 

-Zo 

To  illustrate  mirror  ambiguous  orbits  graphically,  one  periodic  true  relative  orbit  under 
HCW  dynamics  is  chosen  with  initial  condition  shown  in  Table  5.  Figure  14(a)  illustrates 
the  three  mirror  ambiguous  orbits  along  with  the  true  relative  orbit  ( (t):  red  line  ,  (a):  blue 
line,  ( b ):  green  line,  (c):  black  line).  It  is  easy  to  see  that  the  four  orbits  generate  the  same 
range  history  pit)  due  to  the  symmetry  of  the  orbits.  The  y-z  projection  of  the  four  orbits 
shown  in  Figure  14(d)  can  easily  distinguish  different  ambiguous  orbits  by  the  inclination 
of  the  relative  orbit  with  respect  to  the  y- axis  in  the  y-z  plane  (indicated  by  the  slope  of  the 
orbit  in  the  y-z  projection)  and  the  offset  in  the  y  direction,  i.e.,  yd o- 
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Table  5.  Initial  Condition  for  True  Relative  Orbit 


Xo  (km) 

Vo  (km) 

z0  (km) 

x0  (m/s) 

y0  (m/s) 

z0  (m/s) 

-35.51 

12.45 

0.065 

0.66 

—2uxq 

39.43 

x  104 


(c)  x-z  projection 


x  104 


(d)  y-z  projection 


Figure  14.  True  Relative  Orbit  and  Three  Mirror  Ambiguous  Orbits 


6.3  Deformed  Ambiguous  Orbits 

Unlike  a  mirror  ambiguous  orbit,  a  deformed  ambiguous  orbit  implies  that  the  in-plane 
and  out-of-plane  motion  have  different  amplitudes  than  those  of  the  true  relative  orbit,  or 
| a'e\  7^  \ae\  and  \z'm\  ^  \zm\  because  Eqs.  (32-37)  imply  that  a  mismatch  of  the  ambiguous 
orbit’s  in-plane  and  out-of-plane  motion  magnitude  with  those  of  the  true  relative  orbit  must 
occur  simultaneously.  For  a  drifting  relative  orbit,  from  Eq.  (40),  we  see  that  \xd\  =  \x'd\. 
From  Eqs.  (37)  and  (38),  we  have  /30  =  (3'0,  \aeXd\  =  \a'ex'd\  and  \ae\  =  \a'e\.  Then 
considering  Eqs  (32-36),  it  is  easy  to  find  that  \zm\  =  \z'm\  and  \ydo\  =  \y'do\-  In  fact.  merely 
by  \ae\  =  | a'e\  and  zrn  =  \z’m\,  the  existence  of  a  deformed  ambiguous  orbit  for  the  drifting 
case  is  excluded. 
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Without  loss  of  generality,  let  us  assume  a'eyd0  =  aeyd0  from  Eq.  (41).  Furthermore, 
assume  the  ratio  of  in-plane  scale  of  the  deformed  ambiguous  orbit  to  the  true  orbit  is  k, 
namely, 

a'e 
UdO 

Thus,  from  Eqs.  (35)  and  (36)  we  have 

P'o  =  Po  ,  for  P  E  [0,  27 r)  (54) 

Substituting  Eqs.  (52-53)  into  Eqs.  (33)  and  (34),  we  have 

Az'm2  cos(2t/)q)  =  3(k2  —  l)a2e  cos(2/3)  +  Az^  cos(2/3)  (55) 

4:z'm2  sin(2^o)  =  3(k2  -  l)a2  sin(2/3)  +  sin(2/3)  (56) 

Squaring  both  sides  of  Eqs.  (55)  and  (56),  adding  them  together  and  extracting  the  square 

root  results  in 

4 z!^  =  y/9 (k2  -  1  )2aj  +  16^  +  24 (k2  -  l)a^?2n  cos  (2,00  -  2^0)  (57) 

Finally,  substituting  Eqs.  (52-53)  and  (57)  into  Eq.  (32),  we  end  up  with  an  equation  only 
in  terms  of  k,  i.e., 

g 

5 al  +  +  8^o  =  ^k2a2e  +  —y2d0+  V 9(A:2-  1)\4  +  16^n+  24 (A:2-  l)aeVmcos(2/)0-  2-0o)  (58) 

Defining  AT  =  A;2  >  0,  C  =  5a2  +  4z2n  +  8r/^0,  this  can  be  simplified  as 

C  ~  5Ka2e  -  | y2d0  =  y/9(K  -  1  )2a4  +  164  +  -  l)a24  cos(2/30  -  2r/,0)  (59) 

Squaring  both  sides  of  Eq.  (59)  we  obtain 

(c  -  5 Kal  -  I*)  =  9(/i  -  1)X  +  164  +  24(X  -  1)44  c°s(2fti  -  2*)  (60) 

Finally  a  fourth  order  polynomial  can  be  derived  as 

16 Og/i  4  +  [18(4  —  2Aa2z2x  cos(2/30  —  2ip0)  —  10Ca2]  K3 

+  [-9 a4e  +  24a2z2j  cos(2/30  -  2r/i0)  +  80a2^0  -  16**  +  C2]  A'2  -  16Cj/^,Ar  +  64 yAd0 

Note  that  in  deriving  Eq.  (60)  from  Eq.  (59)  the  solution  set  of  K  is  increased  by  “squaring” 
both  sides  of  Eq.  (60).  We  return  to  this  issue  later  when  we  discuss  the  existence  of 
solutions  for  K.  For  now,  it  is  obvious  that  Eq.  (61)  is  a  quartic  equation  and  admits  four 
roots.  However,  not  all  the  roots  are  valid,  because  a  valid  K  for  deformed  ambiguous 
orbits  must  satisfy  the  following  conditions. 


kae 

(52) 

1 

jVdO 

(53) 
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K  >  0 
K  ±  1 


(62) 

(63) 


In  order  to  check  the  solution  set  of  Eq.  (61),  we  test  it  with  a  numerical  example  with 
initial  conditions  illustrated  in  Table  5.  If  we  transform  these  initials  conditions  into  ROEs 
and  substitute  them  into  Eq.  (61),  the  resulting  fourth  order  polynomial  is  given  by 

K 4  -  2.4336/1 3  +  1.629/1 2  -  0.199 K  +  0.004  =  0  (64) 

Eq.  (64)  has  four  real  and  positive  solutions 

"  I<\ 

K2 
Ks 
.  K* 

From  the  solution  set,  we  discuss  the  following  values  separately: 

1.  K2  =  1,  i.e.,  k2  =  ±1  denotes  the  true  relative  orbit  and  corresponding  mirror 
ambiguous  orbits. 

2.  /l3  =  0.125. . . ,  denotes  a  valid  deformed  ambiguous  orbit.  With  K  =  K:i  and 
k  =  ±y/K^,  we  obtain  the  ambiguous  ROEs  by  following  Eqs.  (52-57)  and  (32-  37).  Then 
transforming  these  ambiguous  ROEs  back  to  initial  ambiguous  Cartesian  coordinates,  we 
simulate  the  ambiguous  orbit  using  HCW  dynamics  as  shown  in  Figure  15.  From  Figures 
15(b)  and  15(c),  we  see  that  the  in-plane  motion  amplitude  of  the  ambiguous  orbit  is  de¬ 
creased  compared  to  that  of  the  true  relative  orbit  while  the  out-of-plane  motion  amplitude 
of  the  ambiguous  orbit  is  enlarged  to  compensate  the  decreased  in-plane  amplitude  and  thus 
yield  an  identical  range  history.  From  our  previous  analysis  on  mirror  ambiguous  orbits,  it 
is  seen  that  whenever  one  deformed  ambiguous  orbit  is  found,  there  exist  three  other  orbits 
“mirrored”  to  this  deformed  ambiguous  orbit  generating  the  same  range  history.  In  other 
words,  one  valid  K  ^  1  implies  four  deformed  ambiguous  orbits.  Figure  15  illustrates 
the  true  relative  orbit  and  these  four  deformed  ambiguous  orbits.  (  In  fact,  two  of  these 
deformed  ambiguous  orbits  are  corresponding  to  the  second  case  a'ey'd 0  =  —afy,n)  from  Eq. 

(41).) 

3.  K\  and  K4  are  not  valid  solutions  for  an  ambiguous  orbit  and  are  called  fake  solutions 
in  this  paper.  The  reason  that  these  two  K  values  are  not  valid  ambiguous  solutions,  as  was 
mentioned  earlier,  is  that  the  solution  set  of  K  is  increased  by  “squaring”  both  sides  of  Eq. 
(59).  To  explain  this,  let  us  define 

f=C-  5  Kal  -  ^vl  (66) 

and 

g  =  y/9 (K  -  l)2aj  +  16 +  24(A'  -  1  )a?ezl1  cos(2 /30  -  2^0)  (67) 


1.283681667559196 

1.000000000000000 

0.125488302965344 

0.024439677975557 
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x  104 


(b)  x-y  projection 


x  104 


(c)  x-z  projection 


x  104 


(d)  y-z  projection 


Figure  15.  True  Relative  Orbit  and  Four  Deformed  Ambiguous  Orbits 


Also,  we  define 


p 

1  true 

=  f-g 

(68) 

■Pfake 

=  f  +  g 

(69) 

P 

1  com 

=  K\f-g2) 

(70) 

It  is  easy  to  see  that  PtrUe  =  0  is  equivalent  to  Eq.  (59),  Pfake  =  0  is  the  surplus  or  “fake” 
polynomial  equation  generated  by  the  squaring  operation,  and  Pcom  =  0  is  equivalent  to 
Eq.  (61).  Figure  16  illustrates  these  three  polynomials  with  respect  to  K.  As  expected, 
Pfake  has  two  intersections  with  Pcom  at  the  two  fake  solutions  K L  and  K4,  while  PtrUe  has 
an  additional  two  intersections  with  Pcom  at  /l2  and  I\  :i .  In  other  words,  both  Pfake  =  0 
and  Ptme  =  0  contribute  two  roots  each  to  Pcom  =  0  and  only  the  roots  from  PtrUe  =  0 
are  physically  valid  for  our  problem.  With  these  considerations,  it  is  clear  that  there  are 
at  most  two  real  roots  of  physical  interest  and  one  of  them  is  K  =  1.  Therefore,  we  can 
conclude  that  for  any  closed  relative  orbit,  there  exists  at  most  one  valid  K  ^  1  representing 
deformed  ambiguous  orbits. 
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x  109 


Figure  16.  Three  Polynomials  and  Corresponding  Solutions 


6.4  Existence  of  Deformed  Ambiguous  Orbits 

Recall  that  the  essence  of  this  problem  is  to  decide  the  existence  of  a  valid  K  >  0, 

K  ^  1  in  Eq.  (59).  From  the  previous  analysis,  we  know  that  this  equation  has  two 
solutions  (because  Pcom  =  0  has  four  solutions  in  total  in  the  context  of  complex  solution). 
Now  we  discuss  the  existence  of  this  ambiguous  K .  For  convenience,  define  a  function 
h(K)  as 

h(K)  =  5 Ka2e  +  0  +  ^9(A'-  1)  V+  16Pm+  24(P-  l)aeVmcos(2/30  -  2^0)  (71) 

Therefore,  Eq.  (59)  is  equivalent  to  “  h(K)  =  C  ”.  It  is  obvious  that  h(K )  has  the  following 
properties: 

1.  h(K)  =  C  has  two  solutions  in  the  complex  space. 

2.  h(K )  is  a  continuous  function  for  K  e  (0,  oo). 

3.  K  =  1  is  one  solution  for  h(K)  =  C. 

4.  h(K)  — >  oo  as  K  — >  0. 

5.  h(K)  — >  oo  as  K  — >  oo. 

With  these  five  properties,  we  can  qualitatively  plot  h(K )  as  Figure  17.  From  these  sub¬ 
figures,  if  h(K )  transversely  crosses  C  as  shown  in  Figures  17(a)  and  17(b),  then  we  can 
guarantee  h(K )  =  C  has  one  deformed  ambiguous  solution  at  K  =  K' .  However,  it  is 
possible  that  h(K)  is  tangent  to  C  at  K  =  1  as  shown  in  Figure  17(c)  where  K  =  1  is  a 
double  root  and  thus  no  deformed  ambiguous  K  exists.  Therefore,  we  need  to  discuss  the 
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Figure  17.  Cases  of  Transverse  and  Tangent  Intersections  for  h(K)  and  C 


condition  for  the  tangent  case.  If  h(K)  is  tangent  to  C  at  K  =  1,  it  simply  means  that  the 
derivative  of  h(K)  at  the  tangent  point  equals  to  0,  i.e., 


dh 

dK 


K=i  =  ~  ZVdo  +  \al  cos(2/30  -  2^0)  =  0 


(72) 


Note  that  the  condition  is  independent  of  zm.  To  test  the  condition  for  no  deformed  am¬ 
biguous  orbits  in  Eq.  (72)  ,  we  set 


Vdo  =  y \j\  +  \  -  2,0o)  =  3.551  x  104m  (73) 

in  the  numerical  simulation  while  keeping  the  remaining  ROEs  the  same  as  before  in  Table 
5.  With  these  initial  conditions,  the  quartic  equation  for  K  is 

A'4  -  3.527AT3  +  4.303A'2  -  2.027A'  +  0.250  =  0  (74) 
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Eq.  (74)  has  four  real  and  positive  solutions 


"  Kx  ' 

'  1.340247349090013  ' 

K 2 

1.000000010460630 

K 3 

0.999999989539370 

lu 

0.186546575287586 

It  is  easy  to  see  that  K<i  «  K$  «  1  is  the  repeated  root.  The  same  conclusion  can  also  be 
drawn  from  Figure  18. 


(4)  Pcom 


2  -  ' 


-2 . ' .  .  . . - 
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1 
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1 
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K 

(c)  Pfake 
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K 


(b)  Ptrue 


(d)  h(K) 


Figure  18.  Illustration  of  Tangent  Condition 


6.5  Categorization  of  Ambiguous  Relative  Orbits 

As  discussed  earlier,  we  can  categorize  the  ambiguous  orbits  in  terms  of  initial  Cartesian 
coordinates.  However,  this  criterion  is  not  convenient  when  categorizing  the  various  types 
of  ambiguous  orbits.  Therefore,  this  paper  adopts  the  y-z  plane  projection  to  distinguish 
these  orbits.  Two  parameters  pertaining  to  the  y-z  projection  are  used: 

1.  ydo,  the  offset  in  y  direction. 
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2.  s,  the  slope  of  the  semi-major  axis  of  the  projected  ellipse  in  the  y-z  plane  as  shown 
in  Figure  19.  With  the  non-drifting  solution,  it  can  be  shown  that 

«e  ~  4.  ~  \/<A  ~  2alZm  +  Zm  +  sill2  (/3p  -  ifp) 

2 aezm  sin (/30  -  r/’o) 

The  categorization  of  the  various  types  of  ambiguous  orbits  in  terms  of  ydo  and  s  is  shown 
in  Table  6  with  correlation  between  Cartesian  initial  coordinates  and  y-z  projection.  Also 
note  that  we  require  sign  (yd0)  =  sign  (yd0)  and  sign  (s)  =  sign  (s)  to  decide  the  primary 
deformed  ambiguous  type  (cl),  where  (  )  denotes  deformed  values. 


z 


Figure  19.  Illustration  of  Offset  yd0  and  Slope  s  in  y-z  Projection 


Table  6.  Categorization  of  Ambiguous  Orbits 


Types 

Cartesian  coords. 

y-z  projection 

(t) 

(To,  yo,  z0,  To,  2/o,  Zq) 

VdO 

S 

(a) 

(x0, 2/o,  ~zq,  To,  2/o,  —Zq) 

VdO 

—s 

(b) 

O 

o 

1 

o 

1 

o 

o 

1 

o 

— 2/rfo 

—s 

(c) 

—  (x0,  2/o,  Zo,  T0, 2/o,  zo) 

— 2/rfo 

s 

(d) 

(To,  Vo,  zo,  To,  i/o,  zq) 

VdO 

s 

(e) 

(To,  2/o,  —zo,  To,  Vo,  —zo) 

VdO 

—s 

(/) 

(— T0,  —  Vo,  zo,  —To,  —  2/o,  zo) 

—VdO 

—s 

(9) 

—  (To,  2/o,  ^o,  To,  2/o,  z0) 

—VdO 

s 
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6.6  Numerical  Results  for  Ambiguity  Analysis 

In  the  following  discussion,  we  show  numerical  examples  of  convergence  of  a  sequen¬ 
tial  filter  to  ambiguous  relative  orbits  by  applying  drifting  mirror,  non-drifting  mirror  and 
deformed  ambiguous  initial  conditions  to  an  EKF  with  HCW  dynamics.  For  both  drifting 
and  non-drifting  examples,  the  chief  orbit  is  chosen  as  a  circular  orbit  with  a  radius  of  7100 
km,  and  the  orbital  element  differences  are  shown  in  Table  7. 

The  measurement  covariance  matrix  is  assumed  to  be  R  =  202  m2  and  the  sample  time 
is  At  =  10  seconds.  The  initial  states  guess  for  the  EKF  is  X0  =  (x0,  y0,  z0,  x0,  y0,  Zq)1  = 
(x'0,  y'o,  z'o,  x'0,  y'0,  z'0)f,  where  %  —  (£),  (a), . . . ,  (/)  in  Table  6,  and  the  true  initial  condition 
X(l  =  (xq,  y().  Zq,  To,  yo,  Za)1  is  calculated  from  orbital  element  differences  in  Table  7. 
We  note  that  the  range  measurements  were  generated  from  integrating  nonlinear  two-body 
equations  of  motion,  so  that  the  range  history  does  not  exactly  match  that  generated  from 
HCW  dynamics.  Figure  20  illustrates  the  drifting  estimated  orbits  with  type  (t),  type  (a), 
type  (6)  and  type  (c)  initial  conditions,  in  which  the  red  line  denotes  the  true  relative  orbit 
and  blue  line  denotes  estimated  orbits. 


(a)  True  initial  condition 


(b)  Type  (a)  ambiguous  condition 


x104  xIO4 


(c)  Type  (b)  ambiguous  condition  (d)  Type  (c)  ambiguous  condition 

Figure  20.  Drifting  True  and  Mirror  Ambiguous  Orbits  in  an  EKF  Simulation 
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Figure  20(b)  shows  the  inversion  of  the  z-motion  for  type  (a)  orbit,  Figure  20(c)  shows 
the  inversion  of  the  xy-motion  for  type  ( h )  orbit,  and  Figure  20(d)  shows  the  inversion  of 
the  xyz-motion  for  type  (c)  orbit.  From  these  figures,  it  is  clear  that  the  filter  can  still  arrive 
at  ambiguous  orbits  despite  the  fact  that  their  range  history  is  different  than  that  of  the 
measurements. 

Table  7.  Orbital  Element  Differences  for  Drifting  and  Non-drifting  Scenarios 


5a 

5e 

Si 

5Q 

5u> 

5M0 

Drifting 

5  km 

0.005 

o 

CO 

o 

0° 

0° 

O 

H 

o 

Non-drifting 

0  km 

0.005 

o 

CO 

o 

0° 

0° 

O 

o 

For  the  closed  orbit,  Figure  21  shows  the  results  of  applying  the  type  (c)  mirror  am¬ 
biguous  condition  in  the  EKF.  From  the  numerical  results,  the  filter  reproduces  the  type 
(c)  mirror  ambiguous  orbit  except  with  a  small  disturbance  at  the  beginning  because  of  the 
mismatch  in  measurement  model  (full  nonlinear)  and  filter  dynamic  model  (HCW),  process 
noise  and  measurement  noise.  It  is  worthy  to  mention  that  we  also  found  cases  where  the 
EKF  filter  converged  to  smooth  mirror  ambiguous  orbits  of  type  (a)  and  (6). 


(a)  3-D  plot 


x  104 


x  104 


x  104 


(c)  x-z  projection  (d)  y-z  projection 


Figure  21.  Type  (c)  Mirror  Ambiguous  Orbit  in  an  EKF  Simulation 
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Also,  in  order  to  check  whether  the  filter  will  converge  to  a  deformed  ambiguous  orbit, 
we  applied  the  type  (e)  ambiguous  initial  condition  to  an  EKF.  Figure  22  illustrates  the 
corresponding  results.  From  these  plots,  the  estimated  orbit  does  converge  to  the  “pre¬ 
designed”  deformed  ambiguous  orbit  after  a  certain  initial  oscillation. 


x  104 


xIO4  xIO4 


Figure  22.  Type  (e)  Deformed  Ambiguous  Orbit  in  an  EKF  Simulation 


Some  general  results  on  the  increased  observability  obtained  by  using  higher  order  non- 
linearities  in  the  EKF  resulting  from  the  Taylor  expansion  about  the  chief  orbit  have  been 
shown.  Therefore,  we  conduct  a  comprehensive  set  of  numerical  simulations  to  illustrate 
the  ability  of  quadratic,  cubic  and  full  two-body  nonlinearities  in  the  filter  model  to  avoid 
or  resist  the  occurrence  of  ambiguous  orbits  for  range-only  measurements.  We  consider  the 
effects  of  varying  three  independent  characteristics: 

1.  Nonlinear  order  of  dynamic  model  used  in  the  EKF.  Again  four  different  dynamic 
models  (first ,  second,  third  order  and  full  nonlinear  dynamic  models)  are  applied  in 
the  EKF  as  defined  earlier. 
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2.  Size  of  the  relative  orbit.  Since  the  accuracy  of  describing  relative  motion  using  a 
lower  order  dynamic  model  decreases  as  the  relative  separation  is  enlarged,  three 
configuration  scenarios  are  taken  into  account  as  shown  in  Table  8,  where  pavg  is  the 
average  relative  range  over  the  course  of  one  orbit  period. 

Table  8.  Orbital  element  Differences  for  Three  Simulation  Scenarios 


relative  orbit  size 

Sa 

5D 

5u> 

5e 

Si 

5M0 

Pavg 

small 

0 

0 

0 

0.005 

0.3° 

O 

O 

60  km 

medium 

0 

0 

0 

0.025 

1.5° 

0.5° 

300  km 

large 

0 

0 

0 

0.05 

3° 

1° 

600  km 

3.  Geometric  variation  of  the  initial  guess  given  to  the  EKF.  In  order  to  keep  the  shape 
and  size  of  the  relative  orbit  unchanged  for  one  set  of  variations,  we  alter  the  initial 
guess  for  the  EKF  geometrically  using  ROEs.  To  achieve  this  purpose,  ae  and  zrn 
must  remain  constant  and  yd0,  f30,  and  ip0  are  manipulated  to  fulfill  the  variation 
of  initial  condition.  For  example,  from  Eq.  (51)  and  Eq.  (76),  varying  from  type 
(t)  to  type  (c)  mirror  ambiguous  orbit,  we  need  to  keep  the  phase  difference  (30  — 
ipo  constant  for  slope.  Meanwhile,  to  reverse  the  sign  of  in-plane  and  out-of-plane 
motion,  we  need  to  vary  the  phase  angles  from  [f30,  ip0\  to  \/30  +  tt,  vj(}  +  7r]  and  offset 
in  y  direction  from  yd0  to  —  yd0.  Numerically,  the  initial  guess  in  the  EKF  is  set  as 
[A>>  i’0,  Vdo]  =  [An  V'o,  Vdo]  +  l[n,  tt,  -^Vdo],  where  7  is  the  variation  coefficient  with 
increment  of  0.01.  The  variation  strategy  is  shown  in  Figure  23  with  the  3D  view  and 
y-z  projection. 

The  results  of  this  comprehensive  numerical  investigation  are  shown  and  discussed  in 
the  following  section. 


Y 


(a)  3-D  plot  (b)  y-z  plot 

Figure  23.  Relative  Orbits  Resulting  from  Variation  of  Initial  Conditions 
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7.0  RESULTS  AND  DISCUSSION 


The  first  set  of  simulations  explores  the  variation  from  type  (t)  to  type  (c)  mirror  ambigu¬ 
ous  orbit  for  the  small  relative  orbit  scenario.  The  numerical  results  are  shown  in  Figure 
24,  which  is  called  a  “step  plot”  and  in  which  the  horizontal  lines  indicate  convergence  and 
the  slant  lines  denote  non-convergence.  (  Note  the  filter  is  regarded  as  converged  to  type 
(i)  orbit  if  the  resulting  yd0  and  s  of  the  estimated  orbit  are  within  90%  of  those  of  type 
(i)).  From  this  figure,  it  is  seen  that  the  filter  with  HCW  dynamic  model  diverges  earliest 
at  7  =  .26,  the  second  order  model  diverges  next  at  7  =  .33  and  lastly  the  third  order  and 
full  nonlinear  models  at  7  «  .53.  This  means  that  higher  order  dynamic  models  have  larger 
basins  of  attraction  to  the  true  orbit  compared  with  lower  order  dynamic  models,  which 
converge  more  easily  to  ambiguous  orbits. 


(c) 


(f) 

(b) 

CD  v  ' 

Q_ 

~  (a) 
(d) 
(t) 
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Figure  24.  Convergence  of  EKF  Simulations  for  Small  Relative  Orbit  Scenario 


HCW 

2nd 

3rd 

Full  nonlinear 


Another  interesting  observation  from  Figure  24  is  that  instead  of  only  converging  to  type 
(t)  and  (c)  orbits,  which  are  the  endpoints  of  the  variation,  the  filter  also  converges  to  type 
(d),  (a),  ( b )  and  (/)  ambiguous  orbits.  The  value  of  7  where  the  EKF  switches  from  conver¬ 
gence  to  one  ambiguous  orbit  to  convergence  to  the  "next"  ambiguous  orbit,  increases  for 
higher  order  nonlinear  models.  In  other  words,  the  use  of  higher  order  nonlinear  models 
tend  to  delay  or  resist  the  appearance  of  ambiguous  orbits  in  the  EKF. 

The  step  plots  for  the  medium  and  large  scenarios  are  shown  in  Figure  25.  Figure  25(a) 
and  25(b)  are  very  similar  to  Figure  24  in  terms  of  the  resistance  of  the  higher  order  non¬ 
linear  models  to  ambiguous  orbits  except  that  for  each  order  model,  the  range  of  7  over 
which  the  EKF  converges  to  a  particular  ambiguous  relative  orbit  has  changed.  To  observe 
more  carefully  the  effects  of  large  separation,  Figures  26  and  27  show  the  step  plots  for  the 
HCW  and  full  nonlinear  models  with  medium  and  large  configuration  scenarios. 
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(a)  Medium  scenario  (b)  Large  scenario 

Figure  25.  Convergence  of  EKF  for  Medium  and  Large  Relative  Orbit  Scenarios 


From  Figure  26,  for  the  HCW  dynamic  model,  the  divergence  from  the  true  orbit  be¬ 
comes  faster  when  switching  from  small  to  medium,  and  also  from  medium  to  large  rela¬ 
tive  orbit  scenarios.  This  result  agrees  with  our  expectation  because  the  range  profile  in  the 
EKF  is  generated  by  the  full  nonlinear  dynamic  model  and  the  HCW  model  has  less  ability 
to  track  the  real  motion  and  real  range. 
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Figure  26.  Convergence  of  EKF  for  HCW  Model  with  Three  Relative  Orbit  Scenarios 


Small  scenario 
Medium  scenario 


- Large  scenario 


On  the  contrary,  Figure  27  shows  that  the  filter  diverges  later  when  switching  from 
smaller  to  larger  separation  scenarios  for  the  full  nonlinear  model.  The  reason  is  that  even 
though  the  ambiguous  orbits  perfectly  exist  for  HCW  dynamics,  they  are  not  exactly  suit¬ 
able  for  the  higher  order  or  full  nonlinear  dynamics.  In  other  words,  a  full  nonlinear  model 
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can  never  generate  the  exact  same  range  history  of  the  true  orbit  from  a  different  orbit  for 
which  there  always  exist  errors  in  range.  Therefore,  the  full  nonlinear  model  can  utilize 
this  mismatch  in  range  to  distinguish  and  discard  the  ambiguous  orbit.  The  same  conclu¬ 
sion  could  also  be  drawn  by  performing  the  same  analysis  with  the  second  and  third  order 
models  in  the  EKF.  Results  would  look  similar  to  the  plots  in  Figures  26  and  27. 

In  total,  the  use  of  higher  order  nonlinear  models  in  the  EKF  can  resist  the  tendency  to 
converge  to  an  ambiguous  orbit.  This  resistive  effect  becomes  more  obvious  when  the  rel¬ 
ative  separation  of  the  formation  configuration  is  enlarged.  It  is  noted  that  we  also  studied 
other  directions  of  variation  such  as  (t)— »(/),  (t)— »(a)  and  (t)—>(b)  ,  and  this  conclusion 
holds  for  these  variations  as  well. 
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Figure  27.  Convergence  of  EKF  for  Full  Nonlinear  Model  with  Three  Scenarios 


■Small  scenario 
Medium  scenario 
■  Large  scenario 
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8.0  CONCLUSIONS 


In  this  report,  the  observability  of  relative  orbit  estimation  for  circular  chief  orbits  with 
different  orders  of  nonlinearities  in  the  filter  model  was  studied.  Using  the  Lie  derivative 
method,  analytical  observability  conditions  were  obtained  for  relative  orbit  models  with 
different  orders  of  nonlinearities.  The  benefits  of  using  higher  order  dynamic  models  were 
manifested  both  by  the  Lie  derivative  method  and  by  the  numerical  analysis  using  quantita¬ 
tive  measures  of  observability  obtained  from  the  observability  Gramian,  covariance  matrix, 
and  estimation  error  RMS. 

To  explain  the  reason  for  the  appearance  of  ambiguous  orbits  in  a  Hill-Clohessy-Wiltshire 
based  EKF,  the  ambiguous  relative  orbits  generated  from  range-only  measurements  was  an¬ 
alyzed  both  by  analytical  derivation  and  numerical  simulation.  Two  categories  of  ambigu¬ 
ous  orbits  were  shown  to  exist:  mirror  ambiguous  orbits  and  deformed  ambiguous  orbits, 
and  conclusions  were  drawn  regarding  the  types  and  numbers  of  ambiguous  orbits.  Specifi¬ 
cally,  it  is  shown  that:  1)  For  non-drifting  relative  orbits,  three  mirror  ambiguous  orbits  and 
four  deformed  ambiguous  orbits  exist.  The  deformed  ambiguous  orbits  disappear  when  the 
tangent  condition  is  satisfied.  2)  For  drifting  relative  orbits,  only  three  mirror  ambiguous 
orbits  exist.  These  facts  are  a  consequence  of  the  global  unobservability  of  an  estimation 
strategy  based  on  range-only  measurements  and  HCW  dynamics.  Also  when  used  as  a 
dynamic  model  in  a  filter  the  higher  order  nonlinear  models  of  relative  motion  were  shown 
to  have  the  potential  to  resist  the  tendency  of  an  extended  Kalman  filter  to  converge  to  an 
ambiguous  relative  orbit,  particularly  for  large  separation  relative  orbits. 
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LIST  OF  SYMBOLS,  ABBREVIATIONS,  AND  ACRONYMS 


Hill-Clohessy-Wiltshire 

Local  Vertical  Local  Horizontal 
Extended  Kalman  Filter 
Observability  Index 
Estimation  Condition  Number 
Root  Mean  Square 
Relative  Orbital  Elements 

Symbol  Description 

/x  Gravitational  Parameter  of  the  Earth 

n  Mean  Motion  of  Chief  Orbit 

rc  Orbital  Radius  of  Chief  Orbit 

oj  Angular  Velocity  of  Chief  Orbit 

vc  Velocity  of  Chief  Spacecraft 

x,  y,  z  Cartesian  Position  Coordinates  in  Chief’s  LVLH  Frame 

x,  ij,  z  Cartesian  Velocity  Coordinates  in  Chief’s  LVLH  Frame 

p  Range  Vector  From  Chief  to  Deputy 

X  Deputy’s  Position  and  Velocity 

Yd  Position  of  Deputy  Spacecraft  Relative  to  Center  of  the  Earth 

F,  Vector  filed  of  zth  Order  Dynamic  Model  of  Relative  Motion 

r,  ?th  Order  Nonlinear  Terms 

vrei  Velocity  of  Deputy  Relative  to  the  Chief 

A  In-plane  Bearing  Angle 

0  Out-of-plane  Bearing  Angle 

Wc  Gramian  for  Continuous  System 

Wd  Gramian  for  Discrete  System 

<f>  State  Transition  Matrix 

H  Jacobian  Matrix  of  the  Output  Relation 

a  Singular  Value  of  a  Matrix 

a  Semi-major  Axis  of  an  Orbit 

e  Eccentricity  of  an  Orbit 

i  Inclination  of  an  Orbit 

0  Argument  of  Ascending  Node  of  an  Orbit 

u>  Argument  of  Periapsis  of  an  Orbit 

M0  Initial  Mean  Anomaly  of  an  Orbit 

P0  Initial  Covariance  Matrix 

R  Measurement  Covariance  Matrix 

At  Sampling  Time  Interval 

ep  Estimation  Position  Error 

ep  Estimation  Velocity  Error 

ae  Magnitude  of  In-plane  Motion 

zm  Magnitude  of  Out-of-plane  Motion 

x,i  Offset  in  x  Direction 

yd  Offset  in  y  Direction 

(3  In-plane  Phase  Angle 

0  Out-of-plane  Phase  Angle 

s  Slope  of  the  Semi-major  Axis  of  the  Projected  Ellipse  in  y-z  Plane 

7  Geometric  Variation  Coefficient 


HCW 

LVLH 

EKF 

01 

CN 

RMS 

ROE 
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